import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import pchip_interpolate
import pandas as pd

# 设置长格式
# np.set_printoptions(precision=15)


#输入参数
f_start = 1e6   #起始频点
f_end = 1000e6  #终止频点
f_step = 1e6    #步长
cable_path = r'ParabolaTL.csv'  #线缆文件路径
point_num = 10000  #插值点数
rw = 0.25e-3  # 线缆半径
RL = 150  #负载端负载
R0 = 150  #源端负载
E = 1  #电场强度
theta_p = 50  #入射波theta
psi_p = -70   #入射波phi
theta_E = -30  #电场极化方向



# 频率参数
f = np.arange(f_start, f_end+f_step, f_step)


Lenf = len(f)
w = 2 * np.pi * f
c0 = 299792458
beta = w / c0
mu0 = 4 * np.pi * 1e-7
epsilon0 = 1 / (mu0 * c0 ** 2)
j = 1j  # Python中使用1j表示虚数单位

# 导入线缆路径
Rou = pd.read_csv(cable_path, header=None).values  # 替换为实际文件路径

l = len(Rou)
u = np.linspace(1, l, num=point_num+1)   # 更密集的插值点
N = len(u)

P = (N - 1) / (l - 1)

# 三次样条插值
x = pchip_interpolate(np.arange(1, l + 1), Rou[:, 0], u)
y = pchip_interpolate(np.arange(1, l + 1), Rou[:, 1], u)
z = pchip_interpolate(np.arange(1, l + 1), Rou[:, 2], u)


# # 绘制三维路径
# plt.figure(1)
# ax = plt.axes(projection='3d')
# ax.plot3D(x, y, z, 'royalblue', linewidth=1.5)
# ax.set_xlabel('x/m')
# ax.set_ylabel('y/m')
# ax.set_zlabel('z/m')
# ax.grid(True)



# 平面波参数

theta = np.deg2rad(theta_p)
psi = np.deg2rad(psi_p + 90)
eta = np.deg2rad(theta_E + 90)

ex = np.cos(eta) * np.sin(theta)
ey = np.sin(eta) * np.cos(psi) - np.cos(eta) * np.cos(theta) * np.sin(psi)
ez = np.sin(eta) * np.sin(psi) + np.cos(eta) * np.cos(theta) * np.cos(psi)

betax = beta * np.cos(theta)
betay = beta * np.sin(theta) * np.sin(psi)
betaz = beta * np.sin(theta) * np.cos(psi)

# 计算电场分量
Ex_Re = 2 * E * ex * np.cos(betax[:, np.newaxis] * x) * np.cos(
    betay[:, np.newaxis] * y - betaz[:, np.newaxis] * z)
Ex_Im = 2 * E * ex * np.cos(betax[:, np.newaxis] * x) * np.sin(
    betay[:, np.newaxis] * y - betaz[:, np.newaxis] * z)
Ey_Re = -2 * E * ey * np.sin(betax[:, np.newaxis] * x) * np.sin(
    betay[:, np.newaxis] * y - betaz[:, np.newaxis] * z)
Ey_Im = 2 * E * ey * np.sin(betax[:, np.newaxis] * x) * np.cos(
    betay[:, np.newaxis] * y - betaz[:, np.newaxis] * z)
Ez_Re = -2 * E * ez * np.sin(betax[:, np.newaxis] * x) * np.sin(
    betay[:, np.newaxis] * y - betaz[:, np.newaxis] * z)
Ez_Im = 2 * E * ez * np.sin(betax[:, np.newaxis] * x) * np.cos(
    betay[:, np.newaxis] * y - betaz[:, np.newaxis] * z)


# 端接负载
RL = 150
R0 = 150

# 初始化参数
Len = np.zeros(N - 1)
VsL = np.zeros((Lenf, N - 1), dtype=np.complex_)
VsR = np.zeros((Lenf, N - 1), dtype=np.complex_)
phi11 = np.zeros((Lenf, N - 1), dtype=np.complex_)
phi12 = np.zeros((Lenf, N - 1), dtype=np.complex_)
phi21 = np.zeros((Lenf, N - 1), dtype=np.complex_)
phi22 = np.zeros((Lenf, N - 1), dtype=np.complex_)

for i in range(N - 1):
    dx = x[i + 1] - x[i]
    dy = y[i + 1] - y[i]
    dz = z[i + 1] - z[i]
    Len[i] = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)
    xi_x = dx / Len[i]
    xi_y = dy / Len[i]
    xi_z = dz / Len[i]

    # 计算Sigma和Psi
    Sigma_L = (np.sin(beta * Len[i]) - beta * Len[i]) / (beta ** 2 * np.sin(beta * Len[i]))
    Sigma_R = (np.sin(beta * Len[i]) - beta * Len[i] * np.cos(beta * Len[i])) / (beta ** 2 * np.sin(beta * Len[i]))
    Psi_L = (np.cos(beta * Len[i]) - 1) / (beta * np.sin(beta * Len[i]))
    Psi_R = -Psi_L

    # 计算电场梯度和值
    ax = (Ex_Re[:, i + 1] - Ex_Re[:, i]) / Len[i]
    ay = (Ey_Re[:, i + 1] - Ey_Re[:, i]) / Len[i]
    az = (Ez_Re[:, i + 1] - Ez_Re[:, i]) / Len[i]
    bx = Ex_Re[:, i]
    by = Ey_Re[:, i]
    bz = Ez_Re[:, i]
    cx = (Ex_Im[:, i + 1] - Ex_Im[:, i]) / Len[i]
    cy = (Ey_Im[:, i + 1] - Ey_Im[:, i]) / Len[i]
    cz = (Ez_Im[:, i + 1] - Ez_Im[:, i]) / Len[i]
    dx_ = Ex_Im[:, i]
    dy_ = Ey_Im[:, i]
    dz_ = Ez_Im[:, i]

    # 等效耦合源
    VsL_l = xi_x * ((ax + j * cx) * Sigma_L + (bx + j * dx_) * Psi_L) + \
            xi_y * ((ay + j * cy) * Sigma_L + (by + j * dy_) * Psi_L) + \
            xi_z * ((az + j * cz) * Sigma_L + (bz + j * dz_) * Psi_L)
    VsR_l = xi_x * ((ax + j * cx) * Sigma_R + (bx + j * dx_) * Psi_R) + \
            xi_y * ((ay + j * cy) * Sigma_R + (by + j * dy_) * Psi_R) + \
            xi_z * ((az + j * cz) * Sigma_R + (bz + j * dz_) * Psi_R)
    VsL_v = -(Ex_Re[:, i] + j * Ex_Im[:, i]) * x[i]
    VsR_v = -(Ex_Re[:, i + 1] + j * Ex_Im[:, i + 1]) * x[i + 1]
    VsL[:, i] = VsL_l + VsL_v
    VsR[:, i] = VsR_l + VsR_v

    # 单位长度参数
    hav = (x[i] + x[i + 1]) / 2
    L = (mu0 / (2 * np.pi)) * np.log(2 * hav / rw)
    Zc = c0 * L

    # 传输链参数矩阵
    phi11[:, i] = np.cos(beta * Len[i])
    phi12[:, i] = -j * Zc * np.sin(beta * Len[i])
    phi21[:, i] = -j / Zc * np.sin(beta * Len[i])
    phi22[:, i] = np.cos(beta * Len[i])

# 构建级联矩阵
Phi11 = np.zeros((Lenf, N - 1), dtype=np.complex_)
Phi12 = np.zeros((Lenf, N - 1), dtype=np.complex_)
Phi21 = np.zeros((Lenf, N - 1), dtype=np.complex_)
Phi22 = np.zeros((Lenf, N - 1), dtype=np.complex_)

Phi11[:, 0] = phi11[:, -1]
Phi12[:, 0] = phi12[:, -1]
Phi21[:, 0] = phi21[:, -1]
Phi22[:, 0] = phi22[:, -1]

Cons_V = VsR[:, -1] + Phi11[:, 0] * (-VsL[:, -1] + VsR[:, -2])
Cons_I = Phi21[:, 0] * (-VsL[:, -1] + VsR[:, -2])

for i in range(1, N - 2):
    idx = N - 2 - i
    Phi11[:, i] = Phi11[:, i - 1] * phi11[:, idx] + Phi12[:, i - 1] * phi21[:, idx]
    Phi12[:, i] = Phi11[:, i - 1] * phi12[:, idx] + Phi12[:, i - 1] * phi22[:, idx]
    Phi21[:, i] = Phi21[:, i - 1] * phi11[:, idx] + Phi22[:, i - 1] * phi21[:, idx]
    Phi22[:, i] = Phi21[:, i - 1] * phi12[:, idx] + Phi22[:, i - 1] * phi22[:, idx]
    Cons_V += Phi11[:, i] * (-VsL[:, idx] + VsR[:, idx - 1])
    Cons_I += Phi21[:, i] * (-VsL[:, idx] + VsR[:, idx - 1])

# 最后一级
idx = 0
Phi11[:, -1] = Phi11[:, -2] * phi11[:, idx] + Phi12[:, -2] * phi21[:, idx]
Phi12[:, -1] = Phi11[:, -2] * phi12[:, idx] + Phi12[:, -2] * phi22[:, idx]
Phi21[:, -1] = Phi21[:, -2] * phi11[:, idx] + Phi22[:, -2] * phi21[:, idx]
Phi22[:, -1] = Phi21[:, -2] * phi12[:, idx] + Phi22[:, -2] * phi22[:, idx]
Cons_V += Phi11[:, -1] * (-VsL[:, idx])
Cons_I += Phi21[:, -1] * (-VsL[:, idx])


# 输出参数
denominator = (Phi12[:, -1] - Phi11[:, -1] * R0 + RL * Phi21[:, -1] * R0 - RL * Phi22[:, -1])
I0 = (RL * Cons_I - Cons_V) / denominator  #源端电流
V0 = -R0 * I0                              #源端电压
IL = Cons_I + Phi21[:, -1] * V0 + Phi22[:, -1] * I0  #负载端电流
VL = RL * IL                                         #负载端电压

# print('IO',I0)
# print('V0',V0)
# print('IL',IL)
# print('VL',VL)

# # 转换为dBV
# V0_dBV = 20 * np.log10(np.abs(V0))
# VL_dBV = 20 * np.log10(np.abs(VL))
# I0_dBV = 20 * np.log10(np.abs(I0))
# IL_dBV = 20 * np.log10(np.abs(IL))
#
#
#
# # 读取FEKO数据
# data = pd.read_csv('FEKO_Data_Parabola.csv', header=None, skiprows=2).values
# freq = data[:, 0]
# V0_dBV_MoM = data[:, 1]
# VL_dBV_MoM = data[:, 2]
#
# # 绘图
# plt.rcParams['font.sans-serif'] = ['SimHei']
# plt.rcParams['axes.unicode_minus'] = False
#
# plt.figure(2)
# plt.semilogx(f, V0_dBV, label='V$_{0,Analysis}$')
# plt.semilogx(freq, V0_dBV_MoM, label='V$_{0,MoM}$')
# plt.xlabel('频率/Hz', fontsize=20)
# plt.ylabel('|V$_0$|/dBV', fontsize=20)
# plt.grid(True)
# plt.legend()
#
#
# plt.figure(3)
# plt.semilogx(f, VL_dBV, label='V$_{L,Analysis}$')
# plt.semilogx(freq, VL_dBV_MoM, label='V$_{L,MoM}$')
# plt.xlabel('频率/Hz', fontsize=20)
# plt.ylabel('|V$_L$|/dBV', fontsize=20)
# plt.grid(True)
# plt.legend()
#
# plt.tight_layout()
# plt.show()

